library(rstan)

cat("\nCCES 2016\n\n")

# Set to the local directory of the replication data
setwd("~/Downloads/Replication_Data/")
n_cores <- 8
options(mc.cores = n_cores)

L2016 <- readRDS("Data/CCES-L2016.rds")

##### Stack data to feed to Stan

model_data_2016 <- list(n_items = length(unique(L2016$j)),
                        n_respondents = length(unique(L2016$k)),
                        n_obs = nrow(L2016),
                        k = L2016$k, j = L2016$j, y = L2016$y)

##### Compile to C

model_am <- stan_model(file = "AM_Sigma.stan")

##### Sample from posterior

exclude_pars <- c("beta_raw", "alpha_raw", "zeta_raw")

n_items <- length(unique(L2016$j))
inits <- replicate(n_cores, list(list(zeta_raw = c(-1, 1, -1, 1, -1, runif(n_items-7, -2, 2)))))

posterior_cces_2016 <- sampling(model_am, data = model_data_2016,
                                warmup = 1000, iter = 1250, seed = 1,
                                include = FALSE, pars = exclude_pars,
                                init = inits,
                                refresh = 50, save_warmup = FALSE,
                                control = list(adapt_delta = 0.8,
                                               max_treedepth = 10,
                                               stepsize = 1),
                                open_progress = TRUE, thin = 1, chains = n_cores)

saveRDS(posterior_cces_2016, "Fitted_Models/posterior_cces_2016.rds", compress = "xz")
